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Abstract-Under  wide  sense  stationary  uncorrelated 
scattering  (WSSUS)  conditions,  the  signal  spreading  due  to  a 
random  channel  may  be  described  by  the  scattering  function 
(SF).  In  an  active  acoustic  system,  the  received  signal  is 
modeled  as  the  superposition  of  delayed  and  Doppler  spread 
replicas  of  the  transmitted  waveform.  The  SF  completely 
describes  the  second-order  statistics  of  a  WSSUS  channel  and 
can  be  considered  a  density  function  that  characterizes  the 
average  spread  in  delay  and  Doppler  experienced  by  an  input 
signal  as  it  passes  through  the  channel. 

The  SF  and  its  measurement  will  be  reviewed.  An 
estimator  is  proposed  based  on  a  two-dimensional 
autoregressive  (AR)  model  for  the  scattering  function.  In  order 
to  implement  this  estimator  we  derive  the  minimum  mean 
square  error  estimator  of  the  time-varying  frequency  response 
of  a  linear  channel.  Unlike  conventional  Fourier  methods  the 
AR  approach  does  not  suffer  from  the  usual  convolutional 
smoothing  due  to  the  signal  ambiguity  function.  Simulation 
results  are  given. 

I.  INTRODUCTION 

Transmission  channels  which  spread  the  transmitted 
signal  in  both  time  and  frequency  are  commonly  modeled  as 
random,  time-varying,  space-varying  linear  filters.  Temporal 
spread  is  usually  associated  with  multipath  effects  and 
frequency  spread  with  scatterer  motion.  Under  wide  sense 
stationary  and  uncorrelated  scattering  (WSSUS)  conditions 
the  scattering  function  completely  describes  the  second-order 
statistics  of  the  channel.  It  can  be  viewed  as  a  density 
function  giving  the  average  power  modulation  as  a  function 
of  delay  and  Doppler.  The  SF  is  typically  defined  as  being 
independent  of  the  transmitted  signal.  However,  for  the 
underwater  channel  especially,  this  should  be  considered  true 
only  for  signals  of  similar  bandwidth  and  center  frequency. 

II.  PROBLEM  FORMULATION 

A.  Definitions 

We  first  summarize  the  mathematical  models  that  give 
rise  to  the  scattering  function  estimation  problem.  The 
channel  is  modeled  as  a  stochastic  linear  time-varying  system 
with  impulse  response  h(t,X)  where  h(t,X)  describes  the 

response  of  the  system  at  time  t  to  an  impulse  applied  x 
seconds  prior  [  1  ]  [2] .  Therefore,  if  the  input  to  the  channel  is 
a  signal  s(t),  then  the  output,  x(t)  can  be  written  as 

x(t)=  h(t,x)s(t  -  x)dx .  (l) 
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It  is  assumed  that  the  output  is  the  complex  envelope  and 
therefore  both  s(t)  and  h(t,X )  are  complex.  We  also 
assume  that  h(t,X)  is  zero  mean  for  all  I  and  x,  WSS  in  t 
and.  uncorrelated  in  x.  This  embodies  the  usual  WSSUS 
condition  [3],  Taking  the  Fourier  transform  of  h(t,  X)  with 
respect  to  1  yields  the  spreading  function 

S(<j),x)  =  J"  /?(/,  r)exp(-  jlK(f)i)di  (2) 

which  determines  the  amount  of  spread  in  delay  x  and 
frequency  (f)  that  an  input  signal  undergoes  in  passing  through 
a  time-varying  linear  channel  [1],  Solving  for  h{t,X)  and 
substituting  into  (1)  yields 

-x)exip(j27X</>t)dxd</>.  (3) 

We  see  that  x(t)  is  now  represented  as  the  sum  of  delayed  and 
Doppler  shifted  replicas  of  the  transmitted  signal.  Hence  the 
name  spreading  function  for  S(</>,  z)  is  appropriate. 

Because  h(t,  t)  is  WSS  in  t  and  uncorrelated  in  x  it  can 
be  shown  that  S  {(j),  x)  is  uncorrelated  in  both  (f)  and  x  so  that, 
denoting  the  expectation  operator  by  £(•) , 

x)s{</>',x'))  =  £(s(0, r)|2 ]s{(j)'-(p)8{x'-x)  (4) 

where  the  power,  £,(|5' f)|  ),  for  a  given  Doppler  tf)  and 
delay  x  is  defined  as  the  scattering  function 

P(^,r)  =  £(5(^,r)|2)  (5) 

By  noting  that  the  Fourier  transform  of  h{t,  X)  with  respect 
to  X  yields  the  time-varying  frequency  response  (TVFR) 
H(t,  f ),  which  is  WSS  in  both  1  and  f  we  can  define  the 
channel  ACF  as 

Rh{u,v)=  £[/T(u/  )H{t  +  U,f  +v)]  (6) 

Finally,  it  can  be  shown  that  the  scattering  function  is  related 
to  the  channel  ACF  by  the  two-dimensional  Fourier 
transform 
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P{$,t)=  j"  j"  RH{u,v)e-j2^-m)dudv 


(7) 


9s(u,v)=  As(u,v)2 . 


or  equivalently,  P(</> ,  T)  is  the  power  spectral  density  of 
Hit,  f).  The  channel  ACF,  RH(u,v),  is  also  referred  to  as 
the  two-frequency  correlation  function  [4], 

B.  The  Estimation  Problem 

To  estimate  the  scattering  function  it  is  necessary  to 
either  explicitly  or  implicitly  estimate  the  channel  ACF  due 
to  the  relationship  in  (7).  This  is  difficult  since  the  random 
process  H(t,  f)  is  only  observed  via  the  frequency  domain 
equivalent  of  (1) 

x{t)=  ^H{t,  f)s{f)exp{j27tft)df  (8) 

where 


Using  (10),  one  might  be  inclined  to  use  the  unbiased 
estimate 


Rh  {u,v)=  Ax(u,v)/  As(u,v) 


(12) 


as  was  done  in  [6],  However,  the  correlation  estimate 
becomes  infinite  if  a  (u,v)  =  0.  This  places  severe 

restrictions  on  signal  design  for  realizable  signals.  In  [7]  we 
show  that  the  minimum  mean  square  error  (MMSE)  solution 
for  H( t,f),  is 


H{tn,fm)  = 


S‘  {fm)exy{-  j27tfmtn)x{t„) 


(13) 


//(?,/)=£  h(t,  c)exp(-  ]2n[x)d  x 


for  n= 0,  1,  ...,  N- 1,  m= 0,  1,  ...,  M- 1.  Using  (13)  it  can  be 
shown  that  the  correlation  function  estimate  becomes  [7] 


and  S(f)  is  the  Fourier  transform  of  the  signal  s(t)  (and  not 
the  spreading  function  -  the  number  of  arguments  easily 
distinguish  the  two).  If  the  signal  is  bandlimited  and  f M  , 
is  the  Nyquist  frequency  or  higher,  then  (8)  can  be  expressed 
in  discrete  form  as 


M- 1 

x(tn  )==  A/Xf/  (tn ,  fm  )s  (/,„  )exp(  j27Tfmtn )  (9) 

m—0 

for  n=0,l,...,N-l.  There  are  MN  unknown  samples  of  H(t,f) 
but  only  N  observed  data  samples  so  that  a  consistent 
solution  cannot  be  found.  This  is  known  as  an  overspread 
channel  [5], 

The  theoretical  relationship  between  the  correlation 
function  and  the  transmitted  and  received  signals  can  be 
expressed  in  terms  of  the  time-frequency  (T-F) 
autocorrelation  functions  of  the  signal  and  the  received  time 
series,  As(m,v)  and  Ax(u,v )  [6] 


Rh{u,v) 


E[Ax{u,v)\ 

A,(u,v) 


where  the  T-F  ACF  is  defined  as 


(10) 


As(u,v)=  J”  s*(t  -f)s(t  +  f)e  j2midt 

=  [j*{f-i)s{t  +  x)ej2mfdf. 


(11) 


We  note  that  the  signal  ambiguity  function,  6>v(  u,v)  ,  is  the 
magnitude  squared  of  the  T-F  ACF 


Rh(u,v) 


Ax{u,v)A*{u,v) 

(4(o,o))2 


(14) 


which  is  finite  for  all  signals.  Substituting  the  continuous 
time-frequency  version  of  (13)  into  (8)  yields  an  identity 
after  integration  proving  that  this  is  a  valid  solution  to  the 
estimation  problem.  In  fact,  it  can  also  be  shown  that  (13)  is 
the  solution  of  minimum  norm.  We  note  that  the  MMSE 
estimate  of  the  time -varying  frequency  response  (13)  is 
deterministic  in  the  frequency  direction  (dependent  only  on 
the  transmitted  signal)  and  random  in  the  time  direction.  The 
solution  of  (13)  will  be  used  throughout  the  rest  of  the  paper. 

The  AR  Approach 

We  propose  a  parametric  approach  to  scattering  function 
estimation  based  on  autoregressive  spectral  modeling.  Since 
only  a  few  parameters  must  be  estimated  for  the  AR 
approach,  it  often  can  function  well  when  Fourier  based 
methods  cannot.  Since  simulations  and  any  practical 
implementation  must  be  done  on  a  digital  computer,  a 
discussion  of  sampling  requirements  and  assumptions  is 
appropriate.  As  a  result  of  sampling,  the  scattering  function 
can  only  be  estimated  over  the  Nyquist  band.  Thus,  we  make 
the  very  practical  assumptions  that  the  multipath  (delay) 
spread  is  less  than  L  seconds  and  the  Doppler  spread  is  less 
than  B  Hz.  With  these  assumptions  the  scattering  function 
will  be  estimated  over  the  band 


0  <  X  <  L 

-B/2<(f)<B/2 
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To  prevent  aliasing  in  the  SF,  R^{u,v),  must  be  sampled 
on  a  grid  where  AM  <1  IB  and  Av  <1/L.  We  assume  that 
Am  =  1  IB  and  Av  =  1 1 L  . 

The  scattering  function  is  now  written  using  the  sampled 
form  of  (7) 


"  R[0] 

R[-l] 

R[-p,l 

[  a[0] 

cr2e, 

u  1 

R[l] 

R[0] 

■  R[-  (p,  - 1)] 

a[l] 

= 

0 

_R[p,l 

R[p,  -i] 

R[o] 

L&[^i  ] 

0 

BLkt 


k_  l_ 

~b'T 


-j  M 


B  L 


•  (15) 


where 

a[/]  =  [a[i,0]  a[i,  l]  •••  r?[/,/t2]]r 


If  we  ignore  the  scale  factor  1/BL  and  normalize  the  Doppler 

el  =[l 

o  ••• 

Of  {p2  +l)xl 

and  delay  by  letting  /j  =  (j)  /  B  and  /9  =  z  /  L  and 
r'\k ,  l\  =  R H  ,  j~)  this  becomes 

and 

r[cO] 

Au- 1]  •••  Ac-Pi] 

A/„/2)=  Z  Zr'[*,/]exp(-y2^(/,t-/2/)) 

k=- oo  /=- oo 

R[/]  = 

Ac  l] 

o, 

...  ^ 

1 

—  1/ 2  <  fx  <  1/2,  0  <  /2  <  1 

_Ac  Pi] 

Ac  Pi  -  i]  •••  Ac o] 

which  is  the  usual  definition  of  the  power  spectral  density 
(PSD)  except  for  the  sign  change  (-/2).  To  use  standard  AR 
estimation  techniques  [8]  we  must  account  for  this  sign 
change.  Letting  r[k,l\  =  r'\k  ,  (15)  becomes 


p(f  J2)=  £  fir[k,l]exp{-j2x{flk  +  f2l)) 

k=-°o  l=-° o 

which  is  the  usual  definition  of  the  discrete-time  PSD. 
Therefore,  the  usual  methods  of  2-dimensional  AR  spectral 

2 

estimation  may  be  applied  to  find  the  AR  parameters  ( 7U 

and  «[k,/].  The  spectral  estimator  for  an  A  R  [p  \  ,p2]  quarter 
plane  (QP)  model  is  given  by  [8] 


-i27rU\k  +  fi1)) 

(16) 


The  examples  that  will  be  shown  in  this  paper  use  either 
the  2-D  autocorrelation  method  (ACM)  or  the  2-D 
covariance  method  (CM)  as  defined  in  [8], 


The  Autocorrelation  Method  (ACM) 

The  ACM  requires  an  estimate  of  the  ACF,  samples  of 
which  are  plugged  into  the  Yule- Walker  equations  used  to 
estimate  the  AR  parameters  [8],  The  2-D  Yule-Walker 
equations  are 


(17) 


To  estimate  the  AR  parameters,  we  therefore  need  to 
calculate  the  autocorrelation  function  only  at  the  lags  shown 

in  these  equations  using  r\i,  j\=  r'[k,—l^  =  . 

Assuming  ergodicity,  we  estimate  the  elements  of  the 
autocorrelation  function  defined  in  (6)  using 


n= 0  m= 0 


Lj 

(18) 


The  Covariance  Method  ( CM) 

The  covariance  method,  on  the  other  hand,  requires  an 
estimate  of  H(tf)  and  not  its  ACF.  In  [7]  we  show  that  the 
CM  equations  are 


"  c[o,o] 

C[0,l] 

•  c[o,  px  ] 

r  a[o]  1 

ke, 

c[l,0] 

C[l.l]  •• 

•  C[l,Pl] 

a[l] 

— 

0 

C[Pl,0] 

ck,i]  •• 

•  c  [pi,pA_ 

Ad 

0 

where 


C  [k,i\ 


CTHH[i,:,k,0\ 
C  THH[U:,k,l] 

C  nu  [b  p2  ] 


and  the  colon  denotes  the  entire  range  of  the  index,  i.e. 
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CL[0,:,0,0]  = 

\pHH  [o, 0,0,0]  CHH  [0,1, 0,0]  -■  CHH  [0,  p2  ,0,0]] 


and  each  element  is  calculated  using 


M-lN-l-p2 

CHH\i’  j,k,l\=  ^  Y H(m  —  i,n  —  j)H * \m  —  k, n  + 1\ 

m=p{  n=0 


SIMULATION  RESULTS 

In  2-dimensional  AR  spectral  estimation,  all  causal  AR 
models  are  based  on  a  region  of  support  that  is  either  the 
nonsymmetric  half  plane  (NSHP)  or  the  quarter  plane  (QP). 
In  general  only  the  NSHP  will  yield  the  correct  PSD  if  the 
region  of  support  is  infinite.  However,  it  has  been  observed 
from  simulations  that,  for  sinusoidal  signals  in  noise,  spectral 
estimators  based  on  the  NSHP  perform  poorly,  possibly 
because  the  required  model  order  is  too  high  [8],  All  of  the 
results  presented  herein  utilize  a  2D  quarter  plane  (QP) 
autoregressive  (AR)  model.  Estimates  using  the  ACM  and 
CM  are  compared.  A  comparison  of  results  using  the  NSHP 
and  QP  is  beyond  the  scope  of  this  paper  and  is  an  area  of 
future  work. 

To  demonstrate  the  validity  of  this  approach  the  results 
of  a  number  of  simulations  are  presented.  We  will  assume 

all  data  is  sampled  in  delay  and  Doppler  at  intervals  of  Ar 
and  A^,  respectively.  In  the  simulation  we  define  a  known 

scattering  function,  P,  with  maximum  time  spread  L  and 
maximum  Doppler  spread  B.  We  also  define  a  known 
transmit  waveform  with  time  support  T.  The  samples  of  a 
realization  of  the  spreading  function  are  zero  mean  complex 
Gaussian  variables  with  variance  p(kA^,/AT|  or 

sKlAr)=zkl^^)  where  zkl  ~  CN( °’1)' 

The  received  signal  is  calculated  using  a  discrete  version  of 

(3), 

Bl[ 2Aj  LI A„ 

*(A,)=  y  jr  s(kA^,lAMn\-l\)ej2^lAt 

k=5[ 2At)  1=0 

for  0  <  n  <{T  +  L)/ At  .  Note  that  in  this  calculation 

samples  of  the  transmit  waveform  are  needed  over  the  range 
[ -L ,  T+L ].  If  a  transmitted  signal  is  given  over  an  interval 
from  0  to  T  we  zero-pad  outside  the  interval.  For  a  known 
analytical  expression  such  as  a  Gaussian  pulse,  the  signal  is 
calculated  over  the  entire  range. 

The  first  example  is  for  a  known  AR(1,1)  SF  (defined  by 
(16))  with  time  spread  L= 4  s  and  Doppler  spread  support 
B= 4  Hz.  which  is  interrogated  by  a  Gaussian  probe  pulse  of 
duration  T= 0.357  s.  The  bandwidth  of  a  Gaussian  pulse  is 

W  =  I  /(t^2  )  which  is  1.98  Hz  in  this  case.  The  scattering 
function  is  characterized  by  the  AR  coefficients 
a[0,01=  1.0000,  a[0,l]=0.1 854-0. 5706j,  a[l,0]=0.7000, 


a[l,l]=-0. 1298+0. 3994j.  Fig.  1  shows  examples  of  the 
envelopes  of  the  transmitted  and  received  signals  for  this 
case.  Note  that  the  transmitted  signal  used  for  the  analysis 
has  time  support  over  the  range  [-L,T+L\  and  the  received 
signal  has  time  support  only  over  the  range  [0,L].  All 
contour  plots  are  shown  on  identical  axes  and  contours  are 
given  in  dB.  Fig.  2  shows  the  known  SF  and  the  single  ping 
estimates  for  various  estimators.  The  Fourier  estimate  is 
formed  by  calculating  the  2-D  periodogram  of  the  MMSE 
estimate  of  the  TVFR  or  (13).  This  is  followed  by  AR(1,1) 
estimates  using  both  the  CM  and  the  ACM  estimators. 
Clearly  the  AR  estimators  give  higher  resolution,  more 
accurate  estimates  of  the  scattering  function  for  this  simple 
case. 


(a)  |s(n)| 


Fig.  1 .  (a)  Magnitude  of  transmitted  Gaussian  envelope,  (b) 
Magnitude  of  received  signal  envelope. 
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Fig.  2.  (a)  The  ideal  scattering  function  used  in  the  simulation,  (b) 
The  1-ping  Fourier  estimate,  (c)  The  1-ping  AR(1,1)  CM 
estimate,  (d)  The  1-ping  AR(1,1)  ACM  estimate. 


Scatter  plots  of  AR  parameter  locations  for  50 
realizations  of  the  two  AR(1,1)  estimators  are  shown  in  Fig. 
3.  Solid  lines  on  the  graph  are  drawn  from  the  actual  model 
locations  to  the  average  of  the  50  realizations.  In  almost  all 
cases  the  average  location  of  each  parameter  estimate  is 
biased  toward  the  origin.  The  one  notable  exception  is  for 
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the  a[l,l]  coefficient  using  the  ACM.  The  exact  cause  of 
this  bias  is  a  matter  of  future  investigation.  It  is  also  notable 
that  in  this  case  the  ACM  estimates  of  the  a[  1 .0]  coefficient 
have  significantly  less  scatter  than  the  CM  estimates.  The 
average  scattering  function  estimates  for  these  50  single  ping 
realizations  are  shown  in  Fig.  4. 


Scatter  Plot  of  CM  Coefficient  Locations 


If 


#** 


-0.2  0  0.2 
Real  Part  of  Coefficient 


Scatter  Plot  of  ACM  Coefficient  Locations 


Real  Part  of  Coefficient 

Fig.  3.  Scatter  plot  for  50  realizations  of  1-ping  AR  coefficient 
locations  using  both  the  ACM  and  CM  estimators. 
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Fig.  4.  (a)  The  true  scattering  function,  (b)  The  average  of  50 
Fourier  estimates,  (c)  The  average  of  50  AR(1,1)  CM  estimates, 
(d)  The  average  of  50  AR(1,1)  ACM  estimates. 


Although  we  wish  to  estimate  the  scattering  function 
with  a  single  ping,  the  use  of  multiple  pings  will  improve  the 
accuracy  of  the  estimates  if  the  channel  can  be  considered 
stationary  over  the  time  spanned  by  the  multiple  pings.  Fig.  5 
shows  a  similar  scatter  plot  for  a  case  where  three  pings  are 
used  to  form  the  estimate.  Here,  the  MMSE  estimate  of  the 
TVFR  is  calculated  and  the  corresponding  correlation 
functions  ((18)  into  (17)  for  ACM  or  (13)  into  (19)  for  CM) 
for  AR  estimation  is  formed  for  each  ping.  The  correlation 
functions  are  then  averaged  before  finally  calculating  the  AR 
parameters.  We  see  that  the  variance  of  the  estimates  is 
significantly  reduced  although  the  bias  remains.  Although  it 
is  beyond  the  scope  of  this  paper,  this  indicates  that 


multiping  and/or  recursive  estimation  schemes  may  provide 
robust  estimates  in  environments  where  some  stability  may 
be  assumed  from  ping  to  ping. 


Scatter  Plot  of  CM  Coefficient  Locations 
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Scatter  Plot  of  ACM  Coefficient  Locations 
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Fig.  5.  Scatter  plots  of  50  AR[1,1]  coefficient  locations  for  both 
CM  and  ACM  estimators.  3  pings  are  averaged  to  form  each 
estimate. 


Conclusions 

A  novel  method  of  scattering  function  estimation  based 
on  autoregressive  spectral  modeling  has  been  proposed.  The 
current  implementation  of  this  method  uses  the  MMSE 
estimate  of  the  TVFR  given  a  known  input  waveform  and  the 
received  data.  Preliminary  simulation  results  exhibit  promise 
of  obtaining  high-resolution  estimates  of  the  scattering 
function  from  a  single  ping.  The  results  also  indicate  that  the 
correlation  method  may  be  slightly  more  accurate  on  average 
than  the  autocorrelation  method.  However,  no  claims  of 
optimality  can  be  made  regarding  the  current  estimator. 
Attempts  by  these  authors  to  calculate  the  maximum 
likelihood  estimate  using  the  EM  algorithm  have  failed  due 
to  the  extreme  computational  and  storage  requirements  of  the 
algorithm.  Continuing  research  is  focused  on  improving  this 
technique  using  optimal  estimators  and  waveforms  and  the 
use  of  quarter  plane  versus  nonsymmetric  half  plane 
estimators. 
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